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Abstract 

This work originates from a heart's images tracking which is to generate an ap- 
parent continuous motion, observable through intensity variation from one starting 
image to an ending one both supposed segmented. Given two images po and p±, we 
calculate an evolution process p(t, •) which transports po to p\ by using the optimal 
extended optical flow. In this paper we propose an algorithm based on a fixed point 
formulation and a time-space least squares formulation of the mass conservation 
equation for computing the optimal mass transport problem. The strategy is imple- 
mented in a 2D case and numerical results are presented with a first order Lagrange 
finite element, showing the efficiency of the proposed strategy. 
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1 Introduction 



Modern medical imaging modalities can provide a great amount of informa- 
tion to study the human anatomy and physiological functions in both space 
and time. In cardiac Magnetic Resonance Imaging (MRI) for example, several 
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slices can be acquired to cover the heart in 3D and at a collection of discrete 
time samples over the cardiac cycle. From these partial observations, the chal- 
lenge is to extract the heart's dynamics from these input spatio-temporal data 
throughout the cardiac cycle [12], [T3] . 

Image registration consists in estimating a transformation which insures the 
warping of one reference image onto another target image (supposed to present 
some similarity). Continuous transformations are privileged, the sequence of 
transformations during the estimation process is usually not much considered. 
Most important is the final resulting transformation and not the way one image 
will be transformed to the other. Here, we consider a reasonable registration 
process to continuously map the image intensity functions between two images 
in the context of cardiac motion estimation and modeling. 

The aim of this paper is to present, in the context of extended optical flow, an 
algorithm to compute the optimal time dependent transportation plan without 
using Lagrangian techniques. 

The paper is organized as follows. The introduction is ended, by recalling the 
optimal extended optical flow model (OEOF) . In section 2, the algorithm we 
propose is presented. Its convergence is discussed. In section 3 it is proved that 
solutions obtained with the proposed algorithm are solutions to the optimal 
extended optical flow, that is to say to the time dependent optimal mass 
transportation problem. Section 4 deals with numerical results. A 2D cardiac 
medical image is considered. 

1.1 The OEOF method 

Let us denote by p the intensity function, and by v the velocity of the apparent 
motion of brightness pattern. An image sequence is considered via the gray- 
value map p : Q — (0, 1) x — yM. where Q C M. d is a bounded regular domain, 
the support of images, for d = 1,2,3. If image points move according to the 
velocity field v : Q — > W 1 , then gray values p(t, X(t, x)) are constant along 
motion trajectories X(t, x). One obtains the optical flow equation: 



-pit, X(t, x)) = d t p(t, X(t, x)) + {v\ V xP (t, X{t, x)) ) Rd = 0. (1) 



The assumption that the pixel intensity does not change during the movement 
is in some cases too restrictive. A weakened assumption sometimes called ex- 
tended optical flow, can replace the intensity preservation by a mass preser- 
vation condition which reads: 




(2) 
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The previous equations lead to an ill-posed problem for the unknown (p,v). 
Variational formulations or relaxed minimizing problems for computing jointly 
(p, v) have been first proposed in [4| and after by many other authors. Here 
our concern is somewhat different. Finding (p,v) simultaneously is possible by 
solving the optimal mass transport problem (J3])- (SJ) , developed in [5ll6j . 

Let po and pi be the cardiac images between two times arbitrary fixed to 
zero and one, the mathematical problem reads: find p the gray level function 
defined from Q with values in [0, 1] verifying 

d t p{t, x) + div(t> (t, x)p(t, x)) = 0, in (0, 1) x Q 

(3) 

p(0, x) = p (x); p(l, x) = pi{x) 
The velocity function v is determined in order to minimize the functional: 

inf f 1 [ p(t,x)\\v(t,x)\\ 2 dtdx. (4) 
p> v Jo Jn 

Thus we get an image sequence through the gray- value map p. Let us mention 
[3J, for example, where the optimal mass transportation approach is used in 
images processing. For general properties of optimal transportation, the reader 
is referred to the books by C. Villani [14] and L. Ambrosio et al. [2]. 



2 Algorithm for solving the Optimal Extended Optical Flow 



In what follows, let us specify our hypotheses. 

HI Q is a bounded C 2,a domain satisfying the exterior sphere condition. 
H2 pi G C 1,a (Q) for % — 1,2, and po = p\ on dQ. Moreover there exist two 
constants such that < (3 < pi < (3 in Q. 

Let p° G C^QO, 1] x H) be given by p°(t, x) = (1 - t)p (x) + tp^x). We have 
II^P |lco."([o,i]xH) < C(Po,Pi) and d t p°\ dn = 0. 

For each t G [0, 1], our need for problem (J3J)— (SD is a velocity field vanishing on 
dVL. To do so, the following method is used. 

• Compute 

' -div(p"(t, -)V77) = in 

(5) 

p n (t, -)9 n ?7 = 1 on d£l, 

and set C n {t) = J n d t p n r)dx. 
For each t G [0, 1] compute <p n+1 solution to 
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- div(/>»(t, -)V^ +1 ) = d t p n (t, •), in Q 

(6) 

Set v n+1 = V<p n+1 . 

Compute p n+1 , L 2 -least squares solution to 

d t p n+1 (t,x) + div(v n+1 (t 1 x)p n+1 (t,x)) = 0, in (0,1) x Q 
p" +1 (0,x) = po{x); p n+1 (l,x) = Pl (x). 

For each t E [0, 1], since p n {t, ■), and d t p n {t, •) £ C°' Q (n), theorem 6.14 p. 107 of 
[TT] applies, and there exists a unique <y2 n+1 (t, •) £ C 2,a {VL) solution of problem 
([6]). In problem ([6]) the time is a parameter. As the following regularities with 
respect to time are verified: p n £ C 1,a ; d t p n E C°' a ; C n E C°' a . The classical 
C 2,a (Q) a priori estimates for solutions to elliptic problems allow us to prove 
that (f n+1 is a C 0,a function with respect to time. So we have: 

\\f n+ llc°.«([0,l];CV*(n)) ^ ^(l|C n ||c°.«([0,l]) + II^P n Hc°. Q ([0,l]xn))- 

Consider the extension of ip n+1 by C n outside of the domain Q; still denoted 
by if n+l . Since the right hand side of equation vanishes on dQ, this ex- 
tension is regular, and the function v n+l vanish outside Q and belongs to 
C°' a ([0, 1]; C 1,Q: (R 2 )). 

Define the two flows Xl +1 (s,t,x) E C^fO, 1] x [0, 1] x M 2 ;M 2 ) by 



f t, x) = ±v n+1 (s, Xl +1 (s, t, x)) in (0, 1) 

X^ +1 (t,t,x) =x. 



We have the following 

Lemma 2.1 The L 2 -least squares solution to problem (0) zs ^roen by: 

P ^\t,x) = (i-t f i X t T )) 

r \ i J \ ) p n (t,x) /qn 

P^(^ +1 (l,t,x)) W 

Moreover, i/ < §_ < p n < /3 in [0, 1] x tt, then p n+1 E C^QO, 1] x fi), and 
verifies the same property. 

Proof. We have X™ +1 (1 -s, 1 -i, x) = X" +1 (s, i, x) for every (s, i, z) £ [0, 1] x 
[0, 1] x M 2 (see for example pQ). 

Let us express equation (J7J) along the integral curves of equation (|Sj). The 
L 2 -least squares solution to the ordinary differential equation with initial and 
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final conditions reads 

p n+1 (s,X!l +1 (s,t,x))) = (1 - s)e-So di ^ vn+1 ^ x + +1 ^ x ^ dT p (Xl +1 (0,t,x)) 
+sef* di v(v n+1 (T,xZ +1 {T,t,x))) dT pi (x™ +1 {\, t,x)). 

(10) 

Equation ([6]) gives the following expression for the divergence 

div(v n+1 {s,Xl +1 {s,t,x)) = div(w n+1 (s,X" +1 (l - s, l-t,x)) 

= ^Hp n (s,X^\l-s,l-t,x))). (11) 
as 

The representation formula ([9]) is straightforwardly deduced from ([6]). The 
regularity of the function p n+1 is a consequence of the regularity of the flow 
Xl +l . □ 

Let us now consider the convergence of the algorithm ©-([Tj). 

Theorem 2.2 There exist (p, <p) e C\[0, 1] x H) x C°([0, 1]; C 2 (Q)), L 2 -least 
squares solution, respectively solution to 



d t p(t, x) + div(V<p(i, x)p(t, x)) = 0, in (0, 1) x fi 
p(0,z) = p {x); p{l,x) = px{x) infl 



(12) 



-div(p(t, -)Vy?) = •), in ft 

^ = C(t); V<^ = ondft 



(13) 



C {t) defined by: 



- div(p(t, - )Vr/) = in ft 

p(t,-)d n rj = 1 on<9ft (14) 



Proof. Since ||w ||co-([o,i]) + ||^P || c o, Q([0il]x n) is bounded, 1 1 1| ([o,i] (H)) 
and ||v n+1 ||c* ' Q ([o,i];C 1 ' c "(R 2 )) are uniformly bounded in n. 

From lemma [2TT1 there exists a unique p n+1 , the L 2 -least squares solution of 
(J7|). Let us give an estimate for .D 3 X™ +1 . Starting from 

D 1 X^\s, t, x)) = v n+1 (s, Xl +l (s, t, x)), 
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we deduce (see [T]) 

DsD.X^is, t, x) = D 2 v n+1 (s, Xl +1 (s, t, x))D 3 X^ +1 (s, t, x) 

(!5) 

D 3 Xl +1 (t,t,x) = Id. 
Since D z D 1 Xl +1 (s, t, x) = D 1 D- i Xl +l (s 1 t,x) we get 

D 3 X^ +1 ( S ,t,x) = e -i:^ n+ Hr,XX + \r^)))dr I± 

Thus ||-D3f+ +1 ||co,«([o,i] 2 xiR 2 ) is uniformly bounded in n. 
Since we have p]: 

D 2 Xl +1 (s, t, x) = ( (s, t,x) | D 3 ^ +1 (s, t, x) ) 

we obtain a bound for \\D 2 v n+1 ||c°^([o,i] 2 xR2) independent of n. 

From theorem 12.11 we deduce that ||p n+1 || c i,a([o i] x a) * s uniformly bounded. 
Since the embeddings 

C°' a ([0,l];C 2 ' a (n)) -+C°([0, l];C 2 (n))andC 1 ' Q ([0,l] xfi) ^^([0,1] xfi) 

are relatively compact there is a subsequence of (p n , <p n ) solution to (jSJ)-©, still 
denoted by (p n , <p n ) converging to (p, <p) in ^([0, 1] X O) x C°([0, 1]; C 2 (n)), 
and (p, <p) is the solution of (|T2l-( lT4|) provided the boundary conditions to be 
justified. The condition V<p n |ao = is valid for the approximations <p n (since 
the functions can be extended by C n outside of Q). So the convergence in 
C°([0, 1]; C 2 (Q)) yields the condition for the gradient of limit function. For the 
approximations of function p, the formula given in Lemma \2 . 1 1 combined with 
the regularity result show that the boundary conditions are exactly satisfied. 
These conditions are thus valid for the limit function due to the convergence 
in C 1 . □ 

We will show in the next section that the above least squares solution p is in 
fact a classical solution. 



3 Interpretation of solutions to problem (|T2l-(IT4l) 

In this section it is shown that the solution to problem (|T2|) -( !T4|) is a solution 
to the time dependent optimal mass transportation problem. 
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From one hand, remark that tp solution to problem (fT3|) satisfies: 

1 f 1 

<p-C= Argmin -/ \\d t p + 6xv(pVip)\\ 2 H - Hsl) dt. 

^e J L 2 ((0,l);H 1 (n))4 Jo 



Since the functions (p, <p) are sufficiently regular, we have: 

1 f 1 

^eL 2 ((o,i);/f,;(SJ)niJ 2 (si))' 



1 f 1 

<p - C = Argmin - / \\d t p + div(pV^)||£ 2(n) dt. 

ib£L 2 ((0,l):Hi(n)r\H 2 (Q.))4 J ® 



From an other hand, zero is a bound from below of the functional to be 
minimized with respect to (u,ip): 

= 1/0 \\dtp + div(pV(^ - C))\\l m) dt < 



4 JO \\ u tr ' uiv \r v vr^ ^ ) J\\ L 2 (Q) 

\ Jo ll^t" 



Min U^WdtU + diviuV^Wl^dt. 



{^e^CCO.l)^ 1 ^)), u6L 2 ((0,l);i 2 (n)) 
a t it+div(-?iVV'))eL 2 ((0,i);i 2 (^)) 
9tM+div(?iV?/')=0 
VVIan=0 
V>lan=C 
u(0)=po; «(1)=P1 in f!} 

We deduce that (p, <p), solution to problem (TI2]) -(JT1 |) . satisfies 

1 r 1 

{p,<p)= Argmin -J \\d t u + div(uVip)\\ 2 L 2 iQ) dt. (17) 

{^ei 2 ((0,l);ff 1 (n)), «GL 2 ((0,l);i 2 (^)) 
9 t n+div(-«V^))Gi 2 ((0,l);i 2 (^)) 
dtu+div(u\7ip)=0 

v</>| an =o 

</>lan=C* 
«(0)=po; u(l)=pi in SI} 



Lemma 3.1 Let (p, <p) 6e a solution to problem (1121) - 0141) . T/ien satisfies 
(p,<p)= Argmin jf || div(«VV>)||Ir-i(n) ( 18 

{(9 t u+div(uVi/0=0; V^| S q=0; 

V'lan=C;«(o)=po;«(i)=pi »n ^} 



Proof. This is a simple consequence of d t p = — div(pV<p), and of the regularity 
of div(pV<p) which implies || div(pV<p)||.L 2 (fi) = II div(pV<p)||jj-i(n)- D 

Theorem 3.2 Lei (p, <p) 6e solution to problem (|T2"1) - (IT3|) . i/ie existence of 
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which is given in Theorem \2.^ then it satisfies: 

(p, Vy?) = Argmin / / u\\v || 2 dxdt. (19) 

{dtu+div(uv)=0; u(0)=po;u(l)=pi in f2} ■'O Jfl 

Proof. Choose u regular verifying < /3 < u < /3, and for all t G (0, 1) solve 

inf / u\\v\\ 2 dx. (20) 

{v&L 2 (n)d t u+div(uv)=o} Jn 

Let if = .f/o(fi) be equipped with the following inner product: 

(9, if) = J^u(V6\Vtfj) dx, 

which induces a semi-norm which is equivalent to the if 1 -norm. The Riez's 
theorem claims that for the linear continuous form 

C u (ip) =< -div (uv),i/) >h ; h'=< d t u,ij; > H ;H>, 

there is a unique 9 e H such that 

CM = f u(V6\ Vip ) dx, Vip G H. 
Jn 



Therefore v = V6* and problem ( 1201) is reduced to 



inf 

{tp&H, d t u+div(uX7ip)=0, it>\ d n=C} 



J u\\Vil)\\ 2 dx. (21) 



Since 

problem (1211) reads 



/ u\\Vip\\ 2 dx = || div(uVip)\\ 2 Hh 
Jn 



inf ||div(tiW)|||> (22) 

{V-G-ff, St«+div(uV!/')=0, ^\ an =C} 

or 

inf -||5 t w + div(«V^)||^. (23) 

{ipeH, 9 t «+div(uVi/0=O, V|en=C} 4 

Gathering lemma 13.11 with the previous result proves the theorem. □ 



4 Numerical Approximation of the 2D Optimal Extended Optical 
Flow 



The numerical method is based on a finite element time-space L 2 least squares 
formulation (see [7]) of the linear conservation law ((7j). 
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Define v n+1 as 

v n+l = (1,< +1 X +1 )' 
and for a sufficiently regular function ip denned on Q, set 

~ (dp dip dipV 

and 

2 



Let {ipi ■ ■ ■ ip n} be a basis of a space-time finite element subspace 

Vh = {Pi piecewise regular polynomial functions, with tp(0, •) = tp(l, ■) = 0}, 

for example, a brick Lagrange finite element of order one ([5]). Let 11^ be the 
Lagrange interpolation operator. Let also Wh be the finite element subspace 
of Hq(Q), where the basis functions {ipi ■ ■ -i^m} are the traces at t — of basis 
functions {ipi\f =1 . An approximation of problem (JSJ) is: for a discrete sequence 
of time t compute 

jf (p£(t, •) ( V« +1 - C n (t)) I V^) cfe = f u d tP l(t, -)^ h dx ty h e W h , 

(24) 

and define v n+l = V</^ +1 . The L 2 least squares formulation of problem (J7J) is 
defined in the following way. Consider the functional 



J(c) = \ f (div(^ +1 c) + d t pl + div n h ((l - t)p + tp 

Z JQ 



dx dt. 



This functional is convex and coercive in an appropriate anisotropic Sobolev's 
space [TJ. The minimizer of J is p^ +1 — II^ul — t)po + ^Pi) which is the solution 
to the following problem 

/ div(£ n+1 p™ +1 ) • div(£ n+1 rfx cZt = 

JQ 

I (-d tP n h - ^ n h ((i - 1) Po + t Pl ))) ■ d^(v n+1 ?p h ) dxdt (25) 

J Q 

for all ^ G Vh, where 

JV 

i=i 

Thus an approximation of the solution to problem ([7]) is p^ +1 — II/J (1 — t)p + 



The iterative strategy described in Section 2 is used to compute an approxi- 
mated solution, and to reconstruct the systole to diastole images of a slice of 
a left ventricle. Ten time steps have been used to compute the solution, and 
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10000 degrees of freedom for the time-space least squares finite element. The 
approximated fixed point algorithm converges in about 10 iterations with an 
accuracy of about 10 -7 . In the next figure dj the initial image and the final 
image are presented. In the following figure [21 two intermediate times 1/3 and 




Fig. 1. End of diastole of a left ventricular (a), of systole (b) 



2/3 are shown. 




Fig. 2. Time step 3 and 6 



To summarize, in this work, we present a fixed point algorithm for the compu- 
tation of the time dependent optimal mass transportation problem, allowing 
to handle the images tracking problem. The efficiency of the method has been 
tested with a 2D example. 
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